# sample code for multinomial logistic regression
nall <- 1800
m<-4
run<-0
times<-1000
study<-3
betamean<-matrix(0,(3*(m+1)),1)
sd2mean<-matrix(0,(3*(m+1)),1)
HL_m<-rep(0,times)
AUC_m<-rep(0,times)
test_val<-rep(0,times)

while(run<times)
{	
	
	
# generating data with sample size 1800		
	v1<-rnorm(nall)
	v2<-rnorm(nall)
	v3<-rbinom(nall,1,0.5)
	v4<-rbinom(nall,1,0.5)
	

	 
	xall<-matrix(1,nall,m)
	xall<-cbind(v1,v2,v3,v4)
	
	intcept<-matrix(1,nall,1)
	L1<-2*intcept+0.5*xall[,1]+0.5*xall[,2]+0.5*xall[,3]+0.5*xall[,4]
	L2<-3*intcept+2*xall[,1]+2*xall[,2]+2*xall[,3]+2*xall[,4]
	L3<-1*intcept+1*xall[,1]+1*xall[,2]+1*xall[,3]+1*xall[,4]
	

	p1<-exp(L1)/(1+exp(L1)+exp(L2)+exp(L3))
	p2<-exp(L2)/(1+exp(L1)+exp(L2)+exp(L3))
	p3<-exp(L3)/(1+exp(L1)+exp(L2)+exp(L3))
	p4<-1/(1+exp(L1)+exp(L2)+exp(L3))



	yall<-matrix(0,nall,4)

	for(i in 1:nall)
		yall[i,] <- rmultinom(1,1,c(p1[i],p2[i],p3[i],p4[i]))

	x<-xall[1:900,]
	y<-yall[1:900,]
	xtest<-xall[901:1800,]
	ytest<-yall[901:1800,]


    n<-900
	
	# estimation based on multinomial model
	# add intercept column to old x
	intcept<-matrix(1,n,1)
	x<-cbind(intcept,x)
	zero<-matrix(0,n,(m+1))
	tilde_x<-rbind(cbind(x,zero,zero),cbind(zero,x,zero),cbind(zero,zero,x))
	tilde_y<-matrix(0,3*n,1)
	tilde_y[,1]<-c(y[,1],y[,2],y[,3])
	
	beta0<-matrix(-1,3*(m+1),1)
	beta1<-matrix(0,3*(m+1),1)

	k_cen=0
	epsilon<-10^(-6) 
	while(max(abs(beta1-beta0))>epsilon&k_cen<20)
	{ 
	  beta0<-beta1
	  beta01<-beta0[seq(1,(m+1)),]
	  beta02<-beta0[seq((m+2),(2*m+2)),]
	  beta03<-beta0[seq((2*m+3),(3*m+3)),]
	  p1<-exp(x%*%beta01)/(1+exp(x%*%beta01)+exp(x%*%beta02)+exp(x%*%beta03))
	  p2<-exp(x%*%beta02)/(1+exp(x%*%beta01)+exp(x%*%beta02)+exp(x%*%beta03))
	  p3<-exp(x%*%beta03)/(1+exp(x%*%beta01)+exp(x%*%beta02)+exp(x%*%beta03))
	  p<-rbind(p1,p2,p3)
	  w11<-diag(c(p1*(1-p1)))
	  w22<-diag(c(p2*(1-p2)))
	  w33<-diag(c(p3*(1-p3)))
	  w12<-diag(c(-p1*p2))
	  w13<-diag(c(-p1*p3))
	  w23<-diag(c(-p2*p3))
	  tilde_w<-rbind(cbind(w11,w12,w13),cbind(w12,w22,w23),cbind(w13,w23,w33))
	  d<-t(tilde_x)%*%tilde_w%*%tilde_x
	  e<-t(tilde_x)%*%(tilde_y-p)
	  beta1<-beta0+solve(d+diag(0.000001,3*(m+1)))%*%(e)
	  k_cen=k_cen+1
	  cat("\n\nk_cen=",k_cen)
	}
	if(k_cen<20)
	{
		run<-run+1
		cat("\n\nrun=",run)
		betamean<-betamean+beta1

		hat_beta<-beta1
		hat_beta1<-hat_beta[seq(1,(m+1)),]
		hat_beta2<-hat_beta[seq((m+2),(2*m+2)),]
		hat_beta3<-hat_beta[seq((2*m+3),(3*m+3)),]
		p1<-exp(x%*%hat_beta1)/(1+exp(x%*%hat_beta1)+exp(x%*%hat_beta2)+exp(x%*%hat_beta3))
		p2<-exp(x%*%hat_beta2)/(1+exp(x%*%hat_beta1)+exp(x%*%hat_beta2)+exp(x%*%hat_beta3))
		p3<-exp(x%*%hat_beta3)/(1+exp(x%*%hat_beta1)+exp(x%*%hat_beta2)+exp(x%*%hat_beta3))
		p<-rbind(p1,p2,p3)
		w11<-diag(c(p1*(1-p1)))
		w22<-diag(c(p2*(1-p2)))
		w33<-diag(c(p3*(1-p3)))
		w12<-diag(c(-p1*p2))
		w13<-diag(c(-p1*p3))
		w23<-diag(c(-p2*p3))
		tilde_w<-rbind(cbind(w11,w12,w13),cbind(w12,w22,w23),cbind(w13,w23,w33))
		d<-t(tilde_x)%*%tilde_w%*%tilde_x
		cov_matri_2<-solve(d)
		sd2<-sqrt(diag(cov_matri_2))
		sd2mean<-sd2mean+sd2
		
		# Hosmer-Lemeshow test for multinomial estimation
		ntest<-900
		intcept<-matrix(1,ntest,1)
		xtest<-cbind(intcept,xtest)
		p1<-exp(xtest%*%hat_beta1)/(1+exp(xtest%*%hat_beta1)+exp(xtest%*%hat_beta2)+exp(xtest%*%hat_beta3))
		p2<-exp(xtest%*%hat_beta2)/(1+exp(xtest%*%hat_beta1)+exp(xtest%*%hat_beta2)+exp(xtest%*%hat_beta3))
		p3<-exp(xtest%*%hat_beta3)/(1+exp(xtest%*%hat_beta1)+exp(xtest%*%hat_beta2)+exp(xtest%*%hat_beta3))
		p4<-1-(p1+p2+p3)
		
		psort<-p1+p2+p3
		hltable<-cbind(ytest,p1,p2,p3,p4,psort)
		hltable<-hltable[order(hltable[,9]),]
		
		n_obs1<-matrix(0,10,1)
		n_est1<-matrix(0,10,1)
		n_obs2<-matrix(0,10,1)
		n_est2<-matrix(0,10,1)
		n_obs3<-matrix(0,10,1)
		n_est3<-matrix(0,10,1)
		n_obs4<-matrix(0,10,1)
		n_est4<-matrix(0,10,1)
		sub<-as.integer((ntest)/10)
		
		for(i in 1:9)
		{
			temp<-hltable[(1+(i-1)*sub):(i*sub),1]
			n_obs1[i]<-sum(temp)
			
			temp<-hltable[(1+(i-1)*sub):(i*sub),5]
			n_est1[i]<-sum(temp)
			
			temp<-hltable[(1+(i-1)*sub):(i*sub),2]
			n_obs2[i]<-sum(temp)
			
			temp<-hltable[(1+(i-1)*sub):(i*sub),6]
			n_est2[i]<-sum(temp)
			
			temp<-hltable[(1+(i-1)*sub):(i*sub),3]
			n_obs3[i]<-sum(temp)
			
			temp<-hltable[(1+(i-1)*sub):(i*sub),7]
			n_est3[i]<-sum(temp)
			
			temp<-hltable[(1+(i-1)*sub):(i*sub),4]
			n_obs4[i]<-sum(temp)
			
			temp<-hltable[(1+(i-1)*sub):(i*sub),8]
			n_est4[i]<-sum(temp)
		}
		temp<-hltable[(1+9*sub):(ntest),1]
		n_obs1[10]<-sum(temp)	
		
		temp<-hltable[(1+9*sub):(ntest),5]
		n_est1[10]<-sum(temp)
		
		temp<-hltable[(1+9*sub):(ntest),2]
		n_obs2[10]<-sum(temp)	
		
		temp<-hltable[(1+9*sub):(ntest),6]
		n_est2[10]<-sum(temp)
		
		temp<-hltable[(1+9*sub):(ntest),3]
		n_obs3[10]<-sum(temp)	
		
		temp<-hltable[(1+9*sub):(ntest),7]
		n_est3[10]<-sum(temp)
		
		temp<-hltable[(1+9*sub):(ntest),4]
		n_obs4[10]<-sum(temp)	
		
		temp<-hltable[(1+9*sub):(ntest),8]
		n_est4[10]<-sum(temp)
		
		HL1<-sum((n_obs1[n_est1!=0]-n_est1[n_est1!=0])^2/n_est1[n_est1!=0])		
		HL2<-sum((n_obs2[n_est2!=0]-n_est2[n_est2!=0])^2/n_est2[n_est2!=0])
		HL3<-sum((n_obs3[n_est3!=0]-n_est3[n_est3!=0])^2/n_est3[n_est3!=0])
		HL4<-sum((n_obs4[n_est4!=0]-n_est4[n_est4!=0])^2/n_est4[n_est4!=0])
		
		HL<-HL1+HL2+HL3+HL4
		HL_m[run]<-HL
		
		# AUC score for multinomial estimation
		auctable<-cbind(ytest,p1,p2,p3,p4)
		
		# AUC for y=1,2
		auctable12<-auctable[(auctable[,1]==1|auctable[,2]==1),]
		auctable121<-auctable12[,c(1,5)]
		
		
		tot<-auctable121
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc<-sumrank/(length(pone)*length(pzero))
		auc121<-auc

		
		auctable122<-auctable12[,c(2,6)]
		tot<-auctable122
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc<-sumrank/(length(pone)*length(pzero))
		auc122<-auc
		
		# AUC for y=1,3
		auctable13<-auctable[(auctable[,1]==1|auctable[,3]==1),]
		auctable131<-auctable13[,c(1,5)]
		
		tot<-auctable131
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc<-sumrank/(length(pone)*length(pzero))
		auc131<-auc
		
		auctable133<-auctable13[,c(3,7)]
		
		tot<-auctable133
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc<-sumrank/(length(pone)*length(pzero))
		auc133<-auc
		
		# AUC for y=1,4
		auctable14<-auctable[(auctable[,1]==1|auctable[,4]==1),]
		auctable141<-auctable14[,c(1,5)]
		
		tot<-auctable141
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc<-sumrank/(length(pone)*length(pzero))
		auc141<-auc
		
		auctable144<-auctable14[,c(4,8)]
		
		tot<-auctable144
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc<-sumrank/(length(pone)*length(pzero))
		auc144<-auc
		
		# AUC for y=2,3
		auctable23<-auctable[(auctable[,2]==1|auctable[,3]==1),]
		auctable232<-auctable23[,c(2,6)]
		
		tot<-auctable232
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc<-sumrank/(length(pone)*length(pzero))
		auc232<-auc
		
		
		auctable233<-auctable23[,c(3,7)]
		
		tot<-auctable233
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc<-sumrank/(length(pone)*length(pzero))
		auc233<-auc
		
		# AUC for y=2,4
		auctable24<-auctable[(auctable[,2]==1|auctable[,4]==1),]
		auctable242<-auctable24[,c(2,6)]
		
		tot<-auctable242
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc<-sumrank/(length(pone)*length(pzero))
		auc242<-auc
		
		auctable244<-auctable24[,c(4,8)]
		
		tot<-auctable244
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc<-sumrank/(length(pone)*length(pzero))
		auc244<-auc
		
		# AUC for y=3,4
		auctable34<-auctable[(auctable[,3]==1|auctable[,4]==1),]
		auctable343<-auctable34[,c(3,7)]
		
		tot<-auctable343
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc<-sumrank/(length(pone)*length(pzero))
		auc343<-auc
		
		
		auctable344<-auctable34[,c(4,8)]
		
		tot<-auctable344
		totsort<-tot[order(tot[,2]),]
		onesort<-totsort[totsort[,1]==1,]
		zerosort<-totsort[totsort[,1]==0,]
		pone<-onesort[,2]
		pzero<-zerosort[,2]
		tnp<-rep(0,length(pone))

		for(i in 1: length(pone))
		{
		  count_g<-length(pzero[pzero<pone[i]])
		  count_e<-length(pzero[pzero==pone[i]])
		  tnp[i]<-count_g+0.5*count_e
		}

		sumrank<-sum(tnp)
		auc<-sumrank/(length(pone)*length(pzero))
		auc344<-auc
		
		AUC<-(auc121+auc122+auc131+auc133+auc141+auc144+auc232+auc233+auc242+auc244+auc343+auc344)/12
		AUC_m[run]<-AUC
		
		
		# Testing the assumption of independence of irrelevant alternatives
	
		# get partial inf from full model
		beta_f<-hat_beta[seq((2*m+3),(3*m+3)),]
		cov_f<-cov_matri_2[seq((2*m+3),(3*m+3)),seq((2*m+3),(3*m+3))]
		
		# get partial data
		x3<-x[(y[,3]==1),]
		x4<-x[(y[,4]==1),]
		x34<-rbind(x3,x4)
		y3<-y[(y[,3]==1),]
		y4<-y[(y[,4]==1),]
		y34<-matrix(c(y3[,3],y4[,1]),length(c(y3[,3],y4[,1])),1)
		
		# newton method for restricted model
		beta0<-matrix(-1,(m+1),1)
		beta1<-matrix(0,(m+1),1)


		k_cen=0

		while(max(abs(beta1-beta0))>epsilon&(k_cen<20))
		{ 
			beta0<-beta1
			p<-1/(1+exp(-x34%*%beta0))
			w<-diag(c(p*(1-p)))
			d<-t(x34)%*%w%*%x34
			e<-t(x34)%*%(y34-p)
			beta1<-beta0+solve(d+diag(0.000001,(m+1)))%*%(e)
			k_cen=k_cen+1
			cat("\n\nk_cen=",k_cen)
		}
		if(k_cen<20)
		{
		# covariance matrix for restricted model
		hat_beta_r<-beta1
		p<-1/(1+exp(-x34%*%hat_beta_r))
		w<-diag(c(p*(1-p)))
		d<-t(x34)%*%w%*%x34
		cov_matri_r<-solve(d)
		
		
		# chi-square test statistic for the iia test
		test_val[run]<-t(hat_beta_r-beta_f)%*%solve(cov_matri_r-cov_f)%*%(hat_beta_r-beta_f)
		}else 
			test_val[run]<-1000000
		
			
	}
}

	
	


